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INTRODUCTION 


Current detailed gravimetric geoids are obtained from a combination 
of observed free-air gravity anomalies and geoid undulations obtained from 
satellite derived potential coefficients. A number of such geoids have 
been computed on a 5' grid for the so-called North Atlantic Calibration 
area (Marsh and Vincent, 1973; Marsh and Chang, 1977). These geoids have 
been used in the same area to support GEOS-3 altimeter observations which 
are in turn used to monitor the Gulf Stream. 

An examination of the calculations used in the derivation of these 
geoids as well as associated computer codes, has revealed a few minor 
errors in addition to some inconsistencies in the selection of fundamental 
constants. In view of this and with the availability of new 5' gravity 
data from the DMA gravity library in St. Louis, a revised 5' detailed 
gravimetric geoid for the calibration area has been derived. The method 
of calculation and the derivation of the constants used are explained in 
detail in this report. In addition, the errors of the geoid undulations, 
computed using the rms standard deviations cl the gravity anomalies and 
estimates of the potential coefficient errors, are presented. 


FUNDAMENTALS OF GEO ID COMPUTATION 


Stokes' Integral 

A derivation of the Stokes' Integral yielding the geoid undulation, 

N, is presented in Heiskanen and Moritz (1967) . 

N " N ° * I¥g If Ar s(v) d ° (1) 

where R ■ mean earth radius 

G ■ mean earth gravity 
Ag * surface gravity anomalies 
S(F) * Stokes' function 
do • differential area 

The Stokes' function is given by 

S(f) - sinh m " 6 sln( * /2) + 1 - 5 cos f 

- 3 cos f In £sin(?/2) + sin 2 (y/2) 

where 4* » Cos" 1 Jjsin <t> sin 0 * + cos <J> cos cos (A-A') 

and $,A is the geocentric latitude and longitude respectively of the 
computation point and refer to the position of the variable point 

of integration. 

We also have 


da * cos <p ' dtp ’ dX* 

The term, N , contained in Eq. (1) is called the zero-order undulation 
and is given by 


No 


k£M RAg n 
2 GR " 2 G 


( 2 ) 
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where 

k - Grevitetlonal constant 

AM - difference between assumed and actual earth mass 

Ag ■ global average gravity anomalies ' 

o 

If the actual mass and potential of the earth are not accurately 
known, Eq. (1) yields a geoidal surface consistent with the constants 
used but parallel everywhere to the actual geoid and offset by the constant 
amount N^. In general, because AM is not known and most geoid applications 
require only differences between geoid undulations, N q is set to zero. 

The primary difficulty in performing the integration in Eq. (1) is 
that the desired gravity data is usually not available over the entire 
earth. For this reason the undulation is usually written as three components. 

N - N x + N 2 + N 3 

N, is called the satellite height. It is obtained from the potential 
coefficients and is equivalent to integrating the gravity anomalies, also 
obtained from the potential coefficients, over the entire earth. The second 
term, N 2 , is called the surface height and is obtained by integrating 
Eq. (1) over a finite area using the difference between the observed free- 
air gravity anomalies and those derived from the potential coefficients. 

The third term, would represent the extension of over the remainder 
of the globe, but since sufficient gravity data is lacking is usually 
set equal to zero. This latter action is justified since the decreasing 
wej-jht of the Stokes’ function with increasing V will be applied to areas 
which, although increasing with V , will tend to contain equal amounts of 
positive and negative gravity anomalies, thus yielding a negligible con- 
tribution. The computation of and N 2 is now discussed separately. 
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Satellite Height 


In this section the derivation of the geoidal contribution that is 
a function of the potential coefficients, called the satellite height, 
is discussed. Also presented is the technique for obtaining the 
anomalous gravity field, called the satellite Ag, as a function of the 
potential coefficients. 


A generalized form of Bruns' formula can be written 

T - AW 
N - — 


( 3 ) 


where AW - W - U the difference between the potential of the geold 
o o 

and the potential of the reference ellipsoid. T - W - U where W and U 
are the potentials at the computation point of the geoid and ellipsoid 
respectively, y is the value of gravity at the computation point on the 
ellipsoid. 

If M is the actual mass of the earth, we can write for the gravita- 
tional potential 


W = 


vm oo n n 

l + Z (*) I (C nm cos(mX) + S sin(mX)) P nTn (sin$) 

r n=2 r m=o 


u, 2 r 2 2. 

r— cns <p 


(M 


where 


k 

a 

r 

C S 
nm, nm 


universal gravitational constant 

equatorial radius of reference ellipsoid 

geocentric radius of the point 

unnormalized spherical harmonic coefficients 
of the potential 


P (sinij>) * associated Legendre polynomials 
nm 

to = earth's angular rotation rate 
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* 


am 


A 


A 




The potential at any point on a reference ellipsoid can be defined by 


U 



c 2 ,0 (ref) + 



C, o(* ef ) 

^ I 


- 


P 

nm 


(sin$) 


+ 2 


COS <t> 


(5) 


where is the mass within the reference ellipsoid and C 2 g(ref) and 
g(ref) are reference values of the 2,0 and 4,0 terms of the spherical 
harmonic expansion. Traditionally, only these two terms have been used 
to define the reference ellipsoidal potential. 


N 


Substituting these two potentials into Bruns' formula yields 
kM 


N + 

o 


n n 
m-o 


ry n %2 (f) m-o ( C nm cos(ra ^ + S I™ sln < mX >) P nm< 8in *> 


nm 


where we have assumed M ■ M for the 2,0 and 4,0 terms and 

o 

, jL(M-Mq) _ O'p - .UjbI 
"o ry y 


(6) 


(7) 


where N is again the zero order undulation and C* and S* denote the 
o nm nm 

difference between observed and reference harmonic coefficients. In re- 
ality only 0 ,^ 0 and ^ are differences, all other coefficients being 
equal to the observed values. 

It is obvious from Eq. (7) that tne value of N q will only contribute 

to a constant offset, hence the assumption that M * M and U - W is 

000 

often made and N q , if desired, can be added a posterori. The satellite 
height can now be defined as 


kM 

ry 


z 

n = 2 


(?) 


n 

L 

m=o 


( C nm cos ( mA > + S n m s in (mX ) ) 


P nm (8ln *> 


($> 
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Satellite Gravity 

The same generalized form of Bruns' formula will give the anomalous 
value of gravity between that observed on the geoid and the normal 
gravity on the reference ellipsoid , l.e. 


-6_T 2T 2 AW 

6r " r r 


Ag - — — “ + ~ ~~~ (9 ) 

Using Eq. (3) and noting that N ■ N q + and substituting Eq. (7) yields 


T - Y(N q + Nj) + AW - k l M ~ M o ) + Y N 


( 10 ) 


Substituting this Into the definition for Ag and performing the differentia- 
tion yields 

r z 

IcM - n 


+ 7T 2 (?) 2 (C* cos(mA) + s£ m 8in(mA)) P (sln<>) 

XL n-2 r m-o nm n ® 

+ 7^ n I 2 " ? 2 a b 0 ( c nm C08 < nX ) + S nm 8ln(mX) > P nm (sin«t>) 

~ 2 " ^7? Z, (I)" " (C; m cos(mA) + S* sin (mA ) ) P nm (sin<»>) 

* r n-2 1 m-o nra 

. 2AW 

r (11) 


which can be reduced to the more familiar form 

Ag + kM I (n-1) r av° E (c* cos(mA) + S nm sin(m\)) P nm (sin*> (12) 

0 2 n-2 K r J m-o 

r 

where, by analogy to N q , Ag Q can be called the zero-order gravity anomaly 


A 8, 


■k(M-M„) . 2 AW 

F, + — 


( 13 ) 


which for the same reasons is usually set to zero. 


C 2 ,0 * nd C 4,0 


The unnormalized reference values C 2 Q (ref) and Q (ref) are given by 


C 2 >0 (ref) - -J(n-2) 


C 4.0 <« f > 


-J(n*4) 


( 14 ) 


where J as a function of n Is given by (Heiskanen and Moritz, 1967, p.73) 


2n 


(- 1 ) 


n+1 


3e 


2n 


(2n + 1) (In + 3) 


(1 n + 5n C - A ) 


RE 


(15) 


where 


C, A 

R 

E 

a 

b 

e 


moments of inertia 

principle radius of curvature 

(a 2 - b 2 )* 

semi-major axis 

semi-minor axis 

E/a 


It is also found (ibid) 


C-A 

RE 2 


_2 

45 


me 

qo 


(16) 


where 


m 

J 


2 2 

id a b/kM 
E/b 


and q 0 is given by q when E/u - E/b - e' (ibid p. 71). 


I [(' + Tan' 1 (a') - 1,J 


07) 
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The pest procedure has been to expend Tenge') In a power series 
and develop expressions for C2 g(ref) and Q (ref) using the first two 
terms. There appears to be no reason to continue this because with present 
computer techniques e' t Tan”*(e'), and hence q Q can be determined to any 

desired accuracy. From the above definitions C- -(ref) and C. -(ref) 

* 

can be found exactly and used to determine q and ^ 



C 2,0 - C 2 , 0 (ref) 





(ref) 


( 18 ) 
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Surface Height 


If ehe potential field of the earth, ae given in Eq. (4) were known 
to infinite degree and hence inf initesimal resolution, the geold would 
be given exactly by Eq. (6). Unfortunately this series is known, with 
decreasing accuracy, only as far as degree and order 30. This yields a 
resolution or wavelength of 6° which is not sufficient for the detailed 

work that is desired. As a result, the computation of the geoid must be 

supported by surface observations of gravity anomalies. 

The geoid undulations given by Eq. (8) are equivalent to those that 

would be obtained by integrating in Eq. (1) the satellite gravity 

anomalies given by Eq. (12). However, these satellite gravity anomalies 
are contained implicitly in the observed anomalies and to integrate the 
complete observed anomaly would mean accounting for the satellite field 
twice. Hence the satellite field is subtracted from the observed anomaly 
and only the residual Ip used in Stokes' Integral. This geoid undulation, 
called the surface height is then added to the previous satellite height 
to yield the total geoid height. It should be remembered that the 
contributions to the observed gravity anomalies represented by the 
satellite anomaly is the only contribution that can be effectively inte- 
grated over the entire globe. The residual can be integrated only as far 
as the appropriate data exists. 


The gravity residual then is given by 


V 


^®obs ^®sat 


(19) 


where Ag gat is given by the second term in Eq. (12). 

For the purposes of computation the integral in Eq. (1) is recast as 
a double summation, where 


N 2 ( 4> , X ) 


R 

4ttG 


I 

E 

i--I 


J 

E 

j— J 


A *'iJ 


Sij (4')cos<{> |A4> jAX j 


( 20 ) 
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is the residual 


and N 2 ( $,A ) is the surface height at A, and Ag' 
anomaly and (?) is the Stokes' function at the variable of integration 
position , A'. 

In computing 5' detailed geoids, is usually broken into 3 
components 


N_ - N + N. + N (21) 

i a d c 

where N is the primary contribution and represents the integration of 5* 

a 0-0 
residual anomalies, in the present case, over a square area from 0 to 2 

in each direction away from the computat 4 on point. 

includes the integration of re < 1 anomalies on a 15' grid from 

2° to 3° away from the computation point and includes data on a 

1° grid from 3° to 20° away from the computation point. To find the 

surface height at each 5' interval we add to an average value of 

and N using the values of N, and N adjacent to N since, for the 
c be a 

different resolutions, the 5' positions are not always coincident with the 

centers of the areas used to find N. and ri 

b c. 
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Gravity Data 


Values for observed gravity anomalies have been obtained from the 
Defense Mapping Agency (DMA) gravity library at the Aerospace Center 
In St. Louis. They represent averages over 5' intervals of free-air 
point anomalies where the position of the average is at the center of 
the 5' cell and the cell is identified by the geodetic position of its 
northwest corner. Since the Detailed Geoid Computation (DGC) program 
uses the southwest corner as the cell identifier, the positions are 
corrected to the southwest comer to maintain consistency and avoid 
confusion. 

The geodetic position of each cell is converted to geocentric co- 
ordinates and this position is used to compute satellite height, satellite 
Ag, and the residual gravity anomaly. These values are all assigned to 
the corresponding geodetic position. When the Stokes' integral is 
evaluated, the geodetic positions are again converted to geocentric co- 
ordinates to determine the subtended angle V. The surface heights computed 
are also assigned to the geodetic position, and added to the satellite 
heights tb°i*e to give a geoid undulation at the geodetic position. 

The gravity anomaly is given by 
^®obs ®obs ^ 

where the y used by DMA is the value given by the Geodetic Reference System, 
1967, International Association of Geodesy, Special Publications, No. 3, i.e. 

Y » y e ( 1 + f2sin 2 <J> + f 4 sin^<J) ) (22) 

with 

Y e “ 978031.85 

f 2 - 0.005278895 

f 4 * 0.000023462 

<f> , " geodetic latitude 
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These quantities are computed from the adopted constants in the 
following way. 


/ i _ 1 * _ -2 

ab ' 1 2® 14 


, 2 . 


e m 


) 


( 23 ) 


-f + 5 + 1 f - 26. + 15 2 

-m 2 ~7 f ° ~4 W 


(24) 


4 ~~ 


-1 f 
2 


+ ! f " 


f 25) 


These expressions are obtained from series approximations (Heiskanen 

* * 

and Moritz, 1967, p.76), but unlike the case for C ^ q and q, these 

approximations can be tolerated. This is because all the expressions 

containing y are simply scaled by its value, and minor differences will 

* * 

result in very small percentage errors. However, 'imputing C ^ ^ and ^ 
requires the differencing of two nearly identical - antities that differ 
only in their 6th or 7th decimal places and hence are sensitive to trun- 
cating an infinite series. The results of the geoid calculations are 

& & 

also very sensitive to the values of CL and C, _ that are used and these 
J 2,0 4,0 

numbers should be determined as accurately as possible. 

The constants given in the Geodetic Reference System of 1967 and used 
by DMA to obtain gravity anomalies are not the same constants used in this 
report. To be consistent the following correction must be applied. 


A *n 


ew 


Ag obs + Y 1967 


- Y, 


new 


(26) 
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ERROR ANALYSIS 


The accuracy of the geoid depends on the accuracies of the potential 
coefficients and the gravity data that are used in the geoid computation. 


Potential Coefficients 

The effect of potential coefficient errors on satellite height can 
be seen by examining one element of the double suanation shorn in Eq. (8). 

N 1 (n ’“ ) " 77 ( 7 5 (C nm cos(raX) + S nm sln{nX)) P „ m < #in ♦> 


The differential change in N, (n,m) due to variations in C and S 

I nm i 

can be written 


* 

nm 


(27) 


AN X (n,m) - ^ ( ~ ) P nn (sin *) (AC* na cos(mX) + AS* njn sin (ml)) (28) 

* * 

If the errors in C and S are given in terms of the rros variations a 

c 

and O fi we can obtain cr^ by squaring both sides of (28) and averaging 
over the number of values that are implied in o c and a g at the location 
n,m in the series. 

* (n - (M) (£)" p n m <8ln ♦ ) ( a c (n ' m) cos2 <■*) + % (n,m) sin 2 (mX) ) (29) 

* * 

The cross term AC AS cos (mX ) sin (mX ) has zero average since 
nm nm 

over many measurements of C and S at particular values of n and m, equal 
amounts of positive and negative products are likely, i.e., it is assumed 
that the errors in C and S are uncorrelated. 

The total error to satellite height is found by adding the squares 
of all the individual errors within the summation. 
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• n 
I Z 
n-2 m-0 


2 

°N j (n , m) 


’ (i&j 

'r y ' 


m 

z 

n-2 



n 

Z 

m-0 



(n,m) 


*» 

cos*" 


(«A) 


+ 


o 2 (n,m) sin 2 


(m A) ) P Z (sin*) 


no 


(30) 


It is also necessary to determine the rms error in the satellite gravity, 
to compute the error in surface height arising from errors in the residual 
gravity. This is obtained in an exactly analogous manner to the satellite 
height errors and is given by 


Ag 


- (kM/r 2 ) 2 


<n-ir <f) 


2n 


n 

v 


(o“ (n,m) cos' ( m A > + 


m-0 


2 

o' - (n ,m) 
s 


s in 


(mA)) P‘ (sin*) 
n m 


(31) 


An additional error arises because the potential coefficients above n-N 

r max 

where in the present case N is 30, are not known and cannot be included 

max 

in the series. This so called error of omission has been estimated by 
Rapp (1973, 1975). Like the satellite height and gravity, this error is 
not a strong function of position and hence can be computed separately a 
posteriori and included with the other errors by the addition of small 
constant . 

Computation of the errors given in Eq. (30) and Fq. (31) requires that 

o (n,m) and o (n,m) be known. According to Rapp (1975) the values of o, 

and o are essentially Independent of m and depend only on the degree n. 
s 

He also shows that percentage errors in C and S increase approximately 
linearly with n. Thus, considerable effort may be spared by approximating 
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s 


the potential coefficient errors as a linear function of degree rather 
than using the variance - covariance matrix of the harmonic coefficients. 

Gravity Anomalies 

The error that arises in the surface height due to errors In the 
gravity residuals can be found from Eq. (20) in an analogous way to the 
satellite height. 

In this case, the contribution to N 2 ($,A) from position 
is given by 


N,U,A).. - __R_ A g r s a *.aA 

* 4*G J XJ 1 J 


(32) 


Differential changes are given by 


AN 2 ($, A) ij 


_R_ A (Ag' ) S (Y> A*. A A 

4ttG J 1 J 


(33) 


and the rms error at <p,X due to Ag errors at $',A' is given by 


°N (> ( * ,A) ij 


' 4 * 3 ' 


2 2 


A8 iJ 


S (V) 2 A* 2 AAj 


(34) 


The errors at (J>,A from all that go into the integration are now 

°N (0,X) " A 2 ' £ 'V, S (,,) i3 (/io i 4X 3 )2 (35) 

2 4 * G im _j. j3c .j ij 

The gravity residual is given by Eq. (19) and its rms error is given by 


Ag 


Ag 


+ o 


obs 


Ag 


(36) 
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where o' 
t 


is given by Eq. (31) 

g 


and a 


are the rma errors associated with 


^®obs 


the free air gravity anomalies 


COMPUTATION OF THE REVISED 5' GEOID 


The constants adopted to compute the geoid described In this report 
are those associated with the GEM-8 potential coefficients (Wagner et al, 
1977). The constants required to derive all the quantities used to compute 
the geoid are; 


a - 6378145. m 

kM - 3.986008* 10 14 m 3 sec" 2 


7. 292 115 1467* 10~ 3 sec" 1 


f - 1./298.255 

W 

and the values for mean earth radius and gravity are 

R * 637100C. m 

G - 979800. milligal 

. *» semi-minor axis and the radius to any point on the reference ellipsoid 
at . given by 

b = a (1-f) 
r = a (1-f sin 2 $ ) 

Using these constants, the reference gravity is given by 

* 978030.960508 (1 + 0.00527894364178 sin 2 .* 

+ 0.0000232957478323 sin 4 .* ) 

Since this formula differs slightly from that given by the Geodetic 
Reference System to which the DMA gravity data is referred, the following 
correction is added to the observed free-air gravity anomalies. 


Kb. 


Ag 


obs 


+ 0.89 - . 0428 sin <* + .163 sin i 


The unnormallzed values of ^(ref) and ^(ref) are 

,-3 


C 2 Q (ref) 1.08264312976*10 


C 4,0 (ref) 


2.37097662187*10 


-6 
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Since the GEM-8 values are 

C- n ■ -1.0826249*10" 3 
z y u 

C, 0 - 1.6200262 *10“ 6 

4,0 

* * 

C_ . and C. A . . 

2,0 4,0 are given by 

•k 

C 2,0 - . 18229758- 10 -7 

* 

C 4,0 * -. 75095042 -10" 6 

A list of all the unnormalized GEM-8 coefficients that were used in these 
geoid calculations is given in Appendix I. 

In order to compute the errors in the satellite height arising from 
potential coefficient errors, the following model for the coefficient 
errors was used. 

o,o = 0 N = 2,3,4 

c s 

o o N-4 N>4 

C ’ S = I* 

Hence, the lowest degree terms are most accurately known and the error 
increases linearly, reaching 100% at N * 20. This is an excellent represen- 
tation of the coefficient errors given by Rapp and Rummel (1975) for GEM-6. 
Since the GEM-8 coefficients have been used, which should be an improvement 
over the GEM-6 coefficients, this error estimate may actually be a little 
conservative. 

The errors in surface height were computed from the rms errors in the 
free-air gravity anomalies that accompanied the gravity data from DMA. 
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Using the constants listed in this section and Appendix I and follow- 
ing the method already described, a revised 5' gravimetric geoid was computed 
using the computer facilities of NASA/Wallops Flight Center. This geoid is 
the sum of the satellite and surface heights which are shown as contour maps 
in Figures 1 and 2, respectively. The geoid undulations are presented as a 
contour map in Figure 3. The contour interval for all these figures is 2 
meters. The errors in the geoid have also been computed according to the 
method described and are shown as a contour map in Figure 4. 
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APPENDIX 


The Unnormalized GEM-8 Potential Coefficients 
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